package javavis.jip2d.functions;

import java.util.Arrays;

import javavis.base.ImageType; 
import javavis.base.JIPException;
import javavis.base.parameter.JIPParamBool;
import javavis.base.parameter.JIPParamInt;
import javavis.jip2d.base.FunctionGroup;
import javavis.jip2d.base.JIPFunction;
import javavis.jip2d.base.JIPImage;
import javavis.jip2d.base.bitmaps.JIPImgBitmap;
import javavis.jip2d.base.geometrics.Edge;
import javavis.jip2d.base.geometrics.JIPGeomEdges;
import javavis.jip2d.base.geometrics.JIPGeomPoint;
import javavis.jip2d.base.geometrics.JIPImgGeometric;
import javavis.jip2d.base.geometrics.Point2D;

/**
*FSusan: Implements the Susan method to detect corner points, and edges detection.
*S.M. Smith and J.M. Brady. SUSAN - a new approach to low level image processing. 
*International Journal of Computer Vision, 23(1):45-78, May 1997.
* This implementation is a Java version of the one from the authors
*Input parameters:
*<li>img: Input image
* <li>geomThresh: Geometric threshold: 1850 for corners, 2650 for edges 
* <li>brightThresh: Brightness Threshold
* <li>corners: Detection of corners or edges
* <li>quick: Apply a quick version of Susan
*</ul>
*Output parameters: 
* This function returns a ImageType.POINT image if corners is activated and 
* ImageType.BYTE if not. 
*/
public class FSusan extends JIPFunction {
	private static final long serialVersionUID = -6325253526988864760L;
	
	private static final int LENGTHTHRESH = 2;
	
	public FSusan() {
		super();
		name = "Susan";
		description = "Applies Susan method to detect corners or edges";
		groupFunc = FunctionGroup.Edges;

		JIPParamInt geomThresh = new JIPParamInt("geomThresh", false, true);
		geomThresh.setDefault(1850);
		geomThresh.setDescription("Geometric threshold: 1850 for corners, 2650 for edges");

		JIPParamInt brightThresh = new JIPParamInt("brightThresh", false, true);
		brightThresh.setDefault(13);
		brightThresh.setDescription("Brightness Threshold");

		JIPParamBool corners = new JIPParamBool("corners", false, true);
		corners.setDefault(true);
		corners.setDescription("Detection of corners or edges");

		JIPParamBool quick = new JIPParamBool("quick", false, true);
		quick.setDefault(false);
		quick.setDescription("Apply a quick version of Susan");
		
		addParam(geomThresh);
		addParam(brightThresh);
		addParam(corners);
		addParam(quick);
	}

	public JIPImage processImg(JIPImage img) throws JIPException {
		if (img.getType() == ImageType.BIT || img instanceof JIPImgGeometric) 
			throw new JIPException("Susan can not be applied to this image format");
		
		boolean corners = getParamValueBool("corners");
		boolean quick = getParamValueBool("quick");
		int brightThresh = getParamValueInt("brightThresh");
		int geomThresh = getParamValueInt("geomThresh");
		
		JIPImgBitmap srcbn;
		if (img.getType() == ImageType.COLOR) {
			JIPFunction fctg = new FColorToGray();
			fctg.setParamValue("gray", "BYTE");
			srcbn = (JIPImgBitmap)fctg.processImg(img);
		}
		else srcbn = (JIPImgBitmap)img;
		
		int []lut = setupBrightnessLUT(brightThresh, true);
		
		if (quick) {
			if (corners)
				return susanCornerQuick(srcbn, geomThresh, lut);
			else {
				JIPImgBitmap midImg = (JIPImgBitmap)JIPImage.newImage(img.getWidth(), img.getHeight(), img.getType());
				JIPImgBitmap edges = susanEdgesSmall(srcbn, midImg, geomThresh, lut);
				edgeThin(edges, midImg);
				return extractEdges(midImg, LENGTHTHRESH);
			}
		}
		else
			if (corners)
				return susanCorner(srcbn, geomThresh, lut);
			else {
				JIPImgBitmap midImg = (JIPImgBitmap)JIPImage.newImage(img.getWidth(), img.getHeight(), img.getType());
				JIPImgBitmap edges = susanEdges(srcbn, midImg, geomThresh, lut);
				edgeThin(edges, midImg);
				return extractEdges(midImg, LENGTHTHRESH);
			}
	}

	private boolean is_maximum(double[] pixels, int width, int col, int row) throws JIPException {
		double value = pixels[col + width*row];
		return (value > pixels[(row-3)*width + col-3]
		        && value > pixels[(row-3)*width + col-2]
		        && value > pixels[(row-3)*width + col-1]
  		        && value > pixels[(row-3)*width + col]
		        && value > pixels[(row-3)*width + col+1]
		        && value > pixels[(row-3)*width + col+2]
		        && value > pixels[(row-3)*width + col+3]
                && value > pixels[(row-2)*width + col-3]
		        && value > pixels[(row-2)*width + col-2]
  		        && value > pixels[(row-2)*width + col-1]
		        && value > pixels[(row-2)*width + col]
  		        && value > pixels[(row-2)*width + col+1]
  		        && value > pixels[(row-2)*width + col+2]
  		        && value > pixels[(row-2)*width + col+3]
                && value > pixels[(row-1)*width + col-3]
		        && value > pixels[(row-1)*width + col-2]
  		        && value > pixels[(row-1)*width + col-1]
		        && value > pixels[(row-1)*width + col]
  		        && value > pixels[(row-1)*width + col+1]
  		        && value > pixels[(row-1)*width + col+2]
  		        && value > pixels[(row-1)*width + col+3]
                && value > pixels[(row)*width + col-3]
		        && value > pixels[(row)*width + col-2]
  		        && value > pixels[(row)*width + col-1]
  		        && value >= pixels[(row)*width + col+1]
  		        && value >= pixels[(row)*width + col+2]
  		        && value >= pixels[(row)*width + col+3]
                && value >= pixels[(row+1)*width + col-3]
		        && value >= pixels[(row+1)*width + col-2]
  		        && value >= pixels[(row+1)*width + col-1]
		        && value >= pixels[(row+1)*width + col]
  		        && value >= pixels[(row+1)*width + col+1]
  		        && value >= pixels[(row+1)*width + col+2]
  		        && value >= pixels[(row+1)*width + col+3]
                && value >= pixels[(row+2)*width + col-3]
		        && value >= pixels[(row+2)*width + col-2]
  		        && value >= pixels[(row+2)*width + col-1]
		        && value >= pixels[(row+2)*width + col]
  		        && value >= pixels[(row+2)*width + col+1]
  		        && value >= pixels[(row+2)*width + col+2]
  		        && value >= pixels[(row+2)*width + col+3]
                && value >= pixels[(row+3)*width + col-3]
		        && value >= pixels[(row+3)*width + col-2]
  		        && value >= pixels[(row+3)*width + col-1]
		        && value >= pixels[(row+3)*width + col]
  		        && value >= pixels[(row+3)*width + col+1]
  		        && value >= pixels[(row+3)*width + col+2]
  		        && value >= pixels[(row+3)*width + col+3]);
	}
	
	private int[] setupBrightnessLUT (int thresh, boolean form) {
		int []lut = new int[513];
		double temp;
		
		for (int k=-256; k<257; k++) {
			temp = (double)k/(double)thresh;
			temp *= temp;
			if (form) 
				temp = temp*temp*temp;
			lut[k+256] = (int)(100.0*Math.exp(-temp));
		}
		return lut;
	}
	
	private int ftoi (double a) {
		if (a < 0) 
			return (int)(a-0.5);
		else
			return (int)(a+0.5);
	}
	
	private JIPImgGeometric susanCorner (JIPImgBitmap image, int geomThresh, int[] lut) throws JIPException {
		int w=image.getWidth();
		int h=image.getHeight();
		double []pixels=image.getAllPixels();
		double []r=new double[pixels.length];
		double []cgx=new double[pixels.length];
		double []cgy=new double[pixels.length];
		int n, p, cp, x, y, sq, xx, yy;
		double c;
		
		for (int i=5;i<h-5;i++)
		    for (int j=5;j<w-5;j++) {
		    	n=100;
		        p=(i-3)*w + j - 1;
		        cp=(int)pixels[i*w+j];

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-3; 

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-5;

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-6;

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        if (n<geomThresh){    /* do this test early and often ONLY to save wasted computation */
		        	p+=2;
		        	n+=lut[256+cp-(int)pixels[p++]];
		        	if (n<geomThresh){
		        		n+=lut[256+cp-(int)pixels[p++]];
		        		if (n<geomThresh){
		        			n+=lut[256+cp-(int)pixels[p]];
		        			if (n<geomThresh){
		        				p+=w-6;

		        				n+=lut[256+cp-(int)pixels[p++]];
		        				if (n<geomThresh){
		        					n+=lut[256+cp-(int)pixels[p++]];
		        					if (n<geomThresh){
		        						n+=lut[256+cp-(int)pixels[p++]];
		        						if (n<geomThresh){
		        							n+=lut[256+cp-(int)pixels[p++]];
		        							if (n<geomThresh){
		        								n+=lut[256+cp-(int)pixels[p++]];
		        								if (n<geomThresh){
		        									n+=lut[256+cp-(int)pixels[p++]];
		        									if (n<geomThresh){
		        										n+=lut[256+cp-(int)pixels[p]];
		        										if (n<geomThresh){
		        											p+=w-5;

		        											n+=lut[256+cp-(int)pixels[p++]];
		        											if (n<geomThresh){
		        												n+=lut[256+cp-(int)pixels[p++]];
		        												if (n<geomThresh){
		        													n+=lut[256+cp-(int)pixels[p++]];
		        													if (n<geomThresh){
		        														n+=lut[256+cp-(int)pixels[p++]];
		        														if (n<geomThresh){
		        															n+=lut[256+cp-(int)pixels[p]];
		        															if (n<geomThresh){
		        																p+=w-3;

		        																n+=lut[256+cp-(int)pixels[p++]];
		        																if (n<geomThresh){
		        																	n+=lut[256+cp-(int)pixels[p++]];
		        																	if (n<geomThresh){
		        																		n+=lut[256+cp-(int)pixels[p]];

		        																		if (n<geomThresh) {
		        																			x=0;y=0;
		        																			p=(i-3)*w + j - 1;

		        																			c=lut[256+cp-(int)pixels[p++]]; x-=c; y-=3*c;
																				            c=lut[256+cp-(int)pixels[p++]]; y-=3*c;
																				            c=lut[256+cp-(int)pixels[p]];   x+=c; y-=3*c;
																				            p+=w-3; 
																				    
																				            c=lut[256+cp-(int)pixels[p++]]; x-=2*c; y-=2*c;
																				            c=lut[256+cp-(int)pixels[p++]]; x-=c; y-=2*c;
																				            c=lut[256+cp-(int)pixels[p++]]; y-=2*c;
																				            c=lut[256+cp-(int)pixels[p++]]; x+=c; y-=2*c;
																				            c=lut[256+cp-(int)pixels[p]];   x+=2*c; y-=2*c;
																				            p+=w-5;
																				    
																				            c=lut[256+cp-(int)pixels[p++]]; x-=3*c; y-=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x-=2*c; y-=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x-=c; y-=c;
																				            c=lut[256+cp-(int)pixels[p++]]; y-=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x+=c; y-=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x+=2*c; y-=c;
																				            c=lut[256+cp-(int)pixels[p]];   x+=3*c; y-=c;
																				            p+=w-6;
																		
																				            c=lut[256+cp-(int)pixels[p++]]; x-=3*c;
																				            c=lut[256+cp-(int)pixels[p++]]; x-=2*c;
																				            c=lut[256+cp-(int)pixels[p]];   x-=c;
																				            p+=2;
																				            c=lut[256+cp-(int)pixels[p++]]; x+=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x+=2*c;
																				            c=lut[256+cp-(int)pixels[p]];   x+=3*c;
																				            p+=w-6;
																				    
																				            c=lut[256+cp-(int)pixels[p++]]; x-=3*c; y+=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x-=2*c; y+=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x-=c; y+=c;
																				            c=lut[256+cp-(int)pixels[p++]]; y+=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x+=c; y+=c;
																				            c=lut[256+cp-(int)pixels[p++]]; x+=2*c; y+=c;
																				            c=lut[256+cp-(int)pixels[p]]; x+=3*c; y+=c;
																				            p+=w-5;
																		
																				            c=lut[256+cp-(int)pixels[p++]]; x-=2*c; y+=2*c;
																				            c=lut[256+cp-(int)pixels[p++]]; x-=c; y+=2*c;
																				            c=lut[256+cp-(int)pixels[p++]]; y+=2*c;
																				            c=lut[256+cp-(int)pixels[p++]]; x+=c; y+=2*c;
																				            c=lut[256+cp-(int)pixels[p]]; x+=2*c; y+=2*c;
																				            p+=w-3;
																		
																				            c=lut[256+cp-(int)pixels[p++]]; x-=c; y+=3*c;
																				            c=lut[256+cp-(int)pixels[p++]]; y+=3*c;
																				            c=lut[256+cp-(int)pixels[p]]; x+=c; y+=3*c;
																		
																				            xx=x*x;
																				            yy=y*y;
																				            sq=xx+yy;
																				            if ( sq > ((n*n)/2) )  {
																				            	sq=sq/2;
																				            	if(yy<sq) {
																				            		double divide=y/Math.abs(x);
																				            		sq=Math.abs(x)/x;
																				            		sq=lut[256+cp-(int)pixels[(i+ftoi(divide))*w+j+sq]] +
																				            		   lut[256+cp-(int)pixels[(i+ftoi(2*divide))*w+j+2*sq]] +
																				            		   lut[256+cp-(int)pixels[(i+ftoi(3*divide))*w+j+3*sq]];
																				            	}
																				            	else if(xx<sq) {
																				            		double divide=x/Math.abs(y);
																				            		sq=Math.abs(y)/y;
																				            		sq=lut[256+cp-(int)pixels[(i+sq)*w+j+ftoi(divide)]] +
																				            		   lut[256+cp-(int)pixels[(i+2*sq)*w+j+ftoi(2*divide)]] +
																				            		   lut[256+cp-(int)pixels[(i+3*sq)*w+j+ftoi(3*divide)]];
																				            	}
																				            	if(sq>290){
																				            		r[i*w+j] = geomThresh-n;
																				            		cgx[i*w+j] = (51*x)/n;
																				            		cgy[i*w+j] = (51*y)/n;
																				            	}
																				            }
		        																		}
		        																	}
		        																}
		        															}
		        														}
		        													}
		        												}
		        											}
		        										}
		        									}
		        								}
		        							}
		        						}
		        					}
		        				}
		        			}
		        		}
		        	}
		        }
		    }

	  	JIPGeomPoint resultpoints = (JIPGeomPoint)JIPImage.newImage(w, h, ImageType.POINT);
	  	/* to locate the local maxima */
	  	for (int i=5;i<h-5;i++)
	  		for (int j=5;j<w-5;j++) {
	  			if (r[i*w+j]>0)  {
	  				if (is_maximum(r, w, j, i)) {
	  					resultpoints.addPoint(new Point2D(j, i));
	  				}
	  			}
	  		}
	  	return resultpoints;
	}
	
	private JIPImgGeometric susanCornerQuick (JIPImgBitmap image, int geomThresh, int[] lut) throws JIPException {
		int w=image.getWidth();
		int h=image.getHeight();
		double []pixels=image.getAllPixels();
		double []r=new double[pixels.length];
		int n, p, cp;
		
		for (int i=7;i<h-7;i++)
		    for (int j=7;j<w-7;j++) {
		    	n=100;
		        p=(i-3)*w + j - 1;
		        cp=(int)pixels[i*w+j];

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-3; 

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-5;

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-6;

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        if (n<geomThresh){    /* do this test early and often ONLY to save wasted computation */
		        	p+=2;
		        	n+=lut[256+cp-(int)pixels[p++]];
		        	if (n<geomThresh){
		        		n+=lut[256+cp-(int)pixels[p++]];
		        		if (n<geomThresh){
		        			n+=lut[256+cp-(int)pixels[p]];
		        			if (n<geomThresh){
		        				p+=w-6;

		        				n+=lut[256+cp-(int)pixels[p++]];
		        				if (n<geomThresh){
		        					n+=lut[256+cp-(int)pixels[p++]];
		        					if (n<geomThresh){
		        						n+=lut[256+cp-(int)pixels[p++]];
		        						if (n<geomThresh){
		        							n+=lut[256+cp-(int)pixels[p++]];
		        							if (n<geomThresh){
		        								n+=lut[256+cp-(int)pixels[p++]];
		        								if (n<geomThresh){
		        									n+=lut[256+cp-(int)pixels[p++]];
		        									if (n<geomThresh){
		        										n+=lut[256+cp-(int)pixels[p]];
		        										if (n<geomThresh){
		        											p+=w-5;

		        											n+=lut[256+cp-(int)pixels[p++]];
		        											if (n<geomThresh){
		        												n+=lut[256+cp-(int)pixels[p++]];
		        												if (n<geomThresh){
		        													n+=lut[256+cp-(int)pixels[p++]];
		        													if (n<geomThresh){
		        														n+=lut[256+cp-(int)pixels[p++]];
		        														if (n<geomThresh){
		        															n+=lut[256+cp-(int)pixels[p]];
		        															if (n<geomThresh){
		        																p+=w-3;

		        																n+=lut[256+cp-(int)pixels[p++]];
		        																if (n<geomThresh){
		        																	n+=lut[256+cp-(int)pixels[p++]];
		        																	if (n<geomThresh){
		        																		n+=lut[256+cp-(int)pixels[p]];

		        																		if (n<geomThresh) 
																				            r[i*w+j] = geomThresh-n;
		        																	}
		        																}
		        															}
		        														}
		        													}
		        												}
		        											}
		        										}
		        									}
		        								}
		        							}
		        						}
		        					}
		        				}
		        			}
		        		}
		        	}
		        }
		    }

	  	JIPGeomPoint resultpoints = (JIPGeomPoint)JIPImage.newImage(w, h, ImageType.POINT);
	  	/* to locate the local maxima */
	  	for (int i=5;i<h-5;i++)
	  		for (int j=5;j<w-5;j++) {
	  			if (r[i*w+j]>0)  {
	  				if (is_maximum(r, w, j, i)) {
	  					resultpoints.addPoint(new Point2D(j, i));
	  				}
	  			}
	  		}
	  	return resultpoints;
	}
	
	private void edgeThin (JIPImgBitmap image, JIPImgBitmap midImg) throws JIPException {
		double[] mid = midImg.getAllPixels();
		double[] r = image.getAllPixels();
		int w=image.getWidth();
		int h=image.getHeight();
		
		int   centre, b01, b12, b21, b10,
		      p1, p2, p3, p4,
		      b00, b02, b20, b22,
		      m, n, a=0, b=0, mp, x, y;
		int []l=new int[9];

		for (int i=4; i<h-4; i++)
		    for (int j=4; j<w-4; j++)
		    	if (mid[i*w+j]<8) {
			        centre = (int)r[i*w+j];
			        /* count number of neighbours */
			        mp=(i-1)*w + j-1;
	
			        n = (mid[mp]<8?1:0) +
		        		(mid[mp+1]<8?1:0) +
		        		(mid[mp+2]<8?1:0) +
			        	(mid[mp+w]<8?1:0) +
			        	(mid[mp+w+2]<8?1:0) +
			        	(mid[mp+w+w]<8?1:0) +
			        	(mid[mp+w+w+1]<8?1:0) +
			        	(mid[mp+w+w+2]<8?1:0);
	
			        /* n==0 no neighbours - remove point */
			        if (n==0)
			        	mid[i*w+j]=100;
	
			        /* n==1 - extend line if I can */
			        /* extension is only allowed a few times - the value of mid is used to control this */
			        if (n==1 && mid[i*w+j]<6) {
			          /* find maximum neighbour weighted in direction opposite the
			             neighbour already present. e.g.
			             have: O O O  weight r by 0 2 3
			                   X X O              0 0 4
			                   O O O              0 2 3     */
	
			            l[0]=(int)r[(i-1)*w+j-1]; l[1]=(int)r[(i-1)*w+j]; l[2]=(int)r[(i-1)*w+j+1];
			            l[3]=(int)r[(i  )*w+j-1]; l[4]=0;                 l[5]=(int)r[(i  )*w+j+1];
			            l[6]=(int)r[(i+1)*w+j-1]; l[7]=(int)r[(i+1)*w+j]; l[8]=(int)r[(i+1)*w+j+1];
	
			            if (mid[(i-1)*w+j-1]<8) { 
			            	l[0]=0; l[1]=0; l[3]=0; l[2]*=2; 
			            	l[6]*=2; l[5]*=3; l[7]*=3; l[8]*=4; 
			            }
			            else { 
			            	if (mid[(i-1)*w+j]<8) { 
			            		l[1]=0; l[0]=0; l[2]=0; l[3]*=2; 
			            		l[5]*=2; l[6]*=3; l[8]*=3; l[7]*=4; 
			            	}
			            	else { 
			            		if (mid[(i-1)*w+j+1]<8) { 
			            			l[2]=0; l[1]=0; l[5]=0; l[0]*=2; 
			            			l[8]*=2; l[3]*=3; l[7]*=3; l[6]*=4; 
			            		}
			            		else { 
			            			if (mid[(i)*w+j-1]<8) { 
			            				l[3]=0; l[0]=0; l[6]=0; l[1]*=2; 
			            				l[7]*=2; l[2]*=3; l[8]*=3; l[5]*=4; 
			            			}
			            			else { 
			            				if (mid[(i)*w+j+1]<8) { 
			            					l[5]=0; l[2]=0; l[8]=0; l[1]*=2; 
			            					l[7]*=2; l[0]*=3; l[6]*=3; l[3]*=4; 
			            				}
			            				else { 
			            					if (mid[(i+1)*w+j-1]<8) { 
			            						l[6]=0; l[3]=0; l[7]=0; l[0]*=2; 
			            						l[8]*=2; l[1]*=3; l[5]*=3; l[2]*=4; 
			            					}
			            					else { 
			            						if (mid[(i+1)*w+j]<8) { 
			            							l[7]=0; l[6]=0; l[8]=0; l[3]*=2; 
			            							l[5]*=2; l[0]*=3; l[2]*=3; l[1]*=4; 
			            						}
			            						else { 
			            							if (mid[(i+1)*w+j+1]<8) { 
			            								l[8]=0; l[5]=0; l[7]=0; l[6]*=2; 
			                                            l[2]*=2; l[1]*=3; l[3]*=3; l[0]*=4; 
			                                        } 
			            						}
			            					}
			            				}
			            			}
			            		}
			            	}
			            }
	
				        m=0;     /* find the highest point */
				        for(y=0; y<3; y++)
				            for(x=0; x<3; x++)
				            	if (l[y+y+y+x]>m) { 
				            		m=l[y+y+y+x]; a=y; b=x; 
				            	}
		
				        if (m>0) {
				            if (mid[i*w+j]<4)
				            	mid[(i+a-1)*w+j+b-1] = 4;
				            else
				            	mid[(i+a-1)*w+j+b-1] = mid[i*w+j]+1;
				            if ( (a+a+b) < 3 ) /* need to jump back in image */ {
				            	i+=a-1;
				            	j+=b-2;
				            	if (i<4) i=4;
				            	if (j<4) j=4;
				            }
				        }
			        }
	
			        /* n==2 */
			        if (n==2) {
			        	/* put in a bit here to straighten edges */
			        	b00 = (mid[(i-1)*w+j-1]<8?1:0); /* corners of 3x3 */
			        	b02 = (mid[(i-1)*w+j+1]<8?1:0);
			        	b20 = (mid[(i+1)*w+j-1]<8?1:0);
			        	b22 = (mid[(i+1)*w+j+1]<8?1:0);
			        	if ( ((b00+b02+b20+b22)==2) && ((b00==1 || b22==1) && (b02==1 || b20==1))) {  
			        		/* case: move a point back into line.
			                	e.g. X O X  CAN  become X X X
			                     	 O X O              O O O
			                     	 O O O              O O O    */
			        		if (b00==1)  {
			        			if (b02==1) { x=0; y=-1; }
			        			else     { x=-1; y=0; }
			        		}
			        		else {
			        			if (b02==1) { x=1; y=0; }
			        			else     { x=0; y=1; }
			        		}
			        		if (r[(i+y)*w+j+x]/centre > 0.7) {
			        			if ( ( (x==0) && (mid[(i+(2*y))*w+j]>7) && (mid[(i+(2*y))*w+j-1]>7) && (mid[(i+(2*y))*w+j+1]>7) ) ||
			        				  ( (y==0) && (mid[(i)*w+j+(2*x)]>7) && (mid[(i+1)*w+j+(2*x)]>7) && (mid[(i-1)*w+j+(2*x)]>7) ) ) {
			        			  mid[(i)*w+j]=100;
			        			  mid[(i+y)*w+j+x]=3;  /* no jumping needed */
			        			}
			        		}
			        	}
			        	else {
			        		b01 = (mid[(i-1)*w+j  ]<8?1:0);
			        		b12 = (mid[(i  )*w+j+1]<8?1:0);
			        		b21 = (mid[(i+1)*w+j  ]<8?1:0);
			        		b10 = (mid[(i  )*w+j-1]<8?1:0);
			        		/* right angle ends - not currently used */
			        		if ( ((b01+b12+b21+b10)==2) && ((b10==1 || b12==1) && (b01==1 || b21==1)) &&
			        				((b01==1 && ((mid[(i-2)*w+j-1]<8) || (mid[(i-2)*w+j+1]<8)))|(b10 ==1 &&((mid[(i-1)*w+j-2]<8) || (mid[(i+1)*w+j-2]<8))) ||
			        				(b12==1 && ((mid[(i-1)*w+j+2]<8) || (mid[(i+1)*w+j+2]<8))) || (b21==1 && ((mid[(i+2)*w+j-1]<8) || (mid[(i+2)*w+j+1]<8)))) ) { 
				        	 /* case; clears odd right angles.
			                 e.g.; O O O  becomes O O O
			                       X X O          X O O
			                       O X O          O X O     */
			        			mid[(i)*w+j]=100;
			        			i--;               /* jump back */
			        			j-=2;
			        			if (i<4) i=4;
			        			if (j<4) j=4;
			        		}
			        	}
			        }
	
			        /* n>2 the thinning is done here without breaking connectivity */
			        if (n>2) {
			            b01 = (mid[(i-1)*w+j  ]<8?1:0);
			            b12 = (mid[(i  )*w+j+1]<8?1:0);
			            b21 = (mid[(i+1)*w+j  ]<8?1:0);
			            b10 = (mid[(i  )*w+j-1]<8?1:0);
			            if ((b01+b12+b21+b10)>1) {
			            	b00 = (mid[(i-1)*w+j-1]<8?1:0);
			            	b02 = (mid[(i-1)*w+j+1]<8?1:0);
						    b20 = (mid[(i+1)*w+j-1]<8?1:0);
						    b22 = (mid[(i+1)*w+j+1]<8?1:0);
				            p1 = (b00==1 || b01==1?1:0);
				            p2 = (b02==1 || b12==1?1:0);
				            p3 = (b22==1 || b21==1?1:0);
				            p4 = (b20==1 || b10==1?1:0);
		
				            int aux=0;
				            if (b01==1 && p2==1) aux++; 
				            if (b12==1 && p3==1) aux++;
				            if (b21==1 && p4==1) aux++;
				            if (b10==1 && p1==1) aux++;
				            if( ((p1 + p2 + p3 + p4) - aux) < 2) {
				            	mid[(i)*w+j]=100;
					            i--;
					            j-=2;
					            if (i<4) i=4;
					            if (j<4) j=4;
				            }
			            }
			        }
			    }

	}
	
	private JIPImgBitmap susanEdges (JIPImgBitmap image, JIPImgBitmap midImg, int geomThresh, int[] lut) throws JIPException {
		int w=image.getWidth();
		int h=image.getHeight();
		double []pixels=image.getAllPixels();
		double []r=new double[pixels.length];
		double []mid=new double[pixels.length];
		Arrays.fill(mid, 100);
		int n, p, cp;
		double z, c;
		int m, a, b, x, y, s;
		boolean do_symmetry;
		
		for (int i=3;i<h-3;i++)
			for (int j=3;j<w-3;j++) {
				n=100;
		        p=(i-3)*w + j - 1;
		        cp=(int)pixels[i*w+j];

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-3; 

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-5;

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-6;

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        if (n<geomThresh){    /* do this test early and often ONLY to save wasted computation */
		        	p+=2;
		        	n+=lut[256+cp-(int)pixels[p++]];
		        	if (n<geomThresh){
		        		n+=lut[256+cp-(int)pixels[p++]];
		        		if (n<geomThresh){
		        			n+=lut[256+cp-(int)pixels[p]];
		        			if (n<geomThresh){
		        				p+=w-6;

		        				n+=lut[256+cp-(int)pixels[p++]];
		        				if (n<geomThresh){
		        					n+=lut[256+cp-(int)pixels[p++]];
		        					if (n<geomThresh){
		        						n+=lut[256+cp-(int)pixels[p++]];
		        						if (n<geomThresh){
		        							n+=lut[256+cp-(int)pixels[p++]];
		        							if (n<geomThresh){
		        								n+=lut[256+cp-(int)pixels[p++]];
		        								if (n<geomThresh){
		        									n+=lut[256+cp-(int)pixels[p++]];
		        									if (n<geomThresh){
		        										n+=lut[256+cp-(int)pixels[p]];
		        										if (n<geomThresh){
		        											p+=w-5;

		        											n+=lut[256+cp-(int)pixels[p++]];
		        											if (n<geomThresh){
		        												n+=lut[256+cp-(int)pixels[p++]];
		        												if (n<geomThresh){
		        													n+=lut[256+cp-(int)pixels[p++]];
		        													if (n<geomThresh){
		        														n+=lut[256+cp-(int)pixels[p++]];
		        														if (n<geomThresh){
		        															n+=lut[256+cp-(int)pixels[p]];
		        															if (n<geomThresh){
		        																p+=w-3;

		        																n+=lut[256+cp-(int)pixels[p++]];
		        																if (n<geomThresh){
		        																	n+=lut[256+cp-(int)pixels[p++]];
		        																	if (n<geomThresh){
		        																		n+=lut[256+cp-(int)pixels[p]];

		        																		if (n<geomThresh) 
																				            r[i*w+j] = geomThresh-n;
		        																	}
		        																}
		        															}
		        														}
		        													}
		        												}
		        											}
		        										}
		        									}
		        								}
		        							}
		        						}
		        					}
		        				}
		        			}
		        		}
		        	}
		        }
		    }

	  	for (int i=4;i<h-4;i++)
	  		for (int j=4;j<w-4;j++) {
	  			if (r[i*w+j]>0)  {
	  		        m=(int)r[i*w+j];
	  		        n=geomThresh - m;
			        cp=(int)pixels[i*w+j];

	  		        if (n>600) {
	  			        p=(i-3)*w + j - 1;
		  		        x=0;y=0;
	
		  		        c=lut[256+cp-(int)pixels[p++]]; x-=c; y-=3*c;
			            c=lut[256+cp-(int)pixels[p++]]; y-=3*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=c; y-=3*c;
			            p+=w-3; 
			    
			            c=lut[256+cp-(int)pixels[p++]]; x-=2*c; y-=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; x-=c; y-=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; y-=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c; y-=2*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=2*c; y-=2*c;
			            p+=w-5;
			    
			            c=lut[256+cp-(int)pixels[p++]]; x-=3*c; y-=c;
			            c=lut[256+cp-(int)pixels[p++]]; x-=2*c; y-=c;
			            c=lut[256+cp-(int)pixels[p++]]; x-=c; y-=c;
			            c=lut[256+cp-(int)pixels[p++]]; y-=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c; y-=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=2*c; y-=c;
			            c=lut[256+cp-(int)pixels[p]];   x+=3*c; y-=c;
			            p+=w-6;
	
			            c=lut[256+cp-(int)pixels[p++]]; x-=3*c;
			            c=lut[256+cp-(int)pixels[p++]]; x-=2*c;
			            c=lut[256+cp-(int)pixels[p]];   x-=c;
			            p+=2;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=2*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=3*c;
			            p+=w-6;
			    
			            c=lut[256+cp-(int)pixels[p++]]; x-=3*c; y+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x-=2*c; y+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x-=c; y+=c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c; y+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=2*c; y+=c;
			            c=lut[256+cp-(int)pixels[p]]; x+=3*c; y+=c;
			            p+=w-5;
	
			            c=lut[256+cp-(int)pixels[p++]]; x-=2*c; y+=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; x-=c; y+=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c; y+=2*c;
			            c=lut[256+cp-(int)pixels[p]]; x+=2*c; y+=2*c;
			            p+=w-3;
	
			            c=lut[256+cp-(int)pixels[p++]]; x-=c; y+=3*c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=3*c;
			            c=lut[256+cp-(int)pixels[p]]; x+=c; y+=3*c;
	
		  		        z = Math.sqrt(x*x + y*y);
		  		        if (z > (0.9*n)) /* 0.5 */ {
		  		            do_symmetry=false;
		  		            if (x==0)
		  		            	z=1000000.0;
		  		            else
		  		            	z=(double)y / x;
		  		            if (z < 0) { z=-z; s=-1; }
		  		            else s=1;
		  		            if (z < 0.5) { /* vert_edge */ a=0; b=1; }
		  		            else { if (z > 2.0) { /* hor_edge */ a=1; b=0; }
		  		            else { /* diag_edge */ if (s>0) { a=1; b=1; }
		  		                                   else { a=-1; b=1; }}}
		  		            if ( (m > r[(i+a)*w+j+b]) && (m >= r[(i-a)*w+j-b]) &&
		  		                 (m > r[(i+(2*a))*w+j+(2*b)]) && (m >= r[(i-(2*a))*w+j-(2*b)]) )
		  		            	mid[i*w+j] = 1;
		  		        }
		  		        else
		  		            do_symmetry=true;
		  		    }
		  		    else 
		  		    	do_symmetry=true;
	
	  		        if (do_symmetry) { 
	  			        p=(i-3)*w + j - 1;
		  		        x=0; y=0; s=0;
	
		  		        /*   |      \
		  		             y  -x-  w
		  		             |        \   */

			            c=lut[256+cp-(int)pixels[p++]]; x+=c; y+=9*c; s+=3*c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=9*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=c; y+=9*c; s-=3*c;
			            p+=w-3; 
			    
			            c=lut[256+cp-(int)pixels[p++]]; x+=4*c; y+=4*c; s+=4*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;   y+=4*c; s+=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=4*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;   y+=4*c; s-=2*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=4*c; y+=4*c; s-=4*c;
			            p+=w-5;
			    
			            c=lut[256+cp-(int)pixels[p++]]; x+=9*c; y+=c; s+=3*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=4*c; y+=c; s+=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;   y+=c; s+=c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;   y+=c; s-=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=4*c; y+=c; s-=2*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=9*c; y+=c; s-=3*c;
			            p+=w-6;
	
			            c=lut[256+cp-(int)pixels[p++]]; x+=9*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=4*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=c;
			            p+=2;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=4*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=9*c;
			            p+=w-6;
			    
			            c=lut[256+cp-(int)pixels[p++]]; x+=9*c; y+=c; s-=3*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=4*c; y+=c; s-=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;   y+=c; s-=c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;   y+=c; s+=c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=4*c; y+=c; s+=2*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=9*c; y+=c; s+=3*c;
			            p+=w-5;
	
			            c=lut[256+cp-(int)pixels[p++]]; x+=4*c; y+=4*c; s-=4*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;   y+=4*c; s-=2*c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=4*c;
			            c=lut[256+cp-(int)pixels[p++]]; x+=c;   y+=4*c; s+=2*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=4*c; y+=4*c; s+=4*c;
			            p+=w-3;
	
			            c=lut[256+cp-(int)pixels[p++]]; x+=c; y+=9*c; s-=3*c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=9*c;
			            c=lut[256+cp-(int)pixels[p]];   x+=c; y+=9*c; s+=3*c;
		  		        
		  		        if (y==0)
		  		        	z = 1000000.0;
		  		        else
		  		        	z = (double)x / y;
		  		        if (z < 0.5) { /* vertical */ a=0; b=1; }
		  		        else { 
		  		        	if (z > 2.0) { /* horizontal */ a=1; b=0; }
		  		        	else { /* diagonal */ 
		  		        		if (s>0) { a=-1; b=1; }
		  		        		else { a=1; b=1; }
		  		        	}
		  		        }
		  		        if ( (m > r[(i+a)*w+j+b]) && (m >= r[(i-a)*w+j-b]) &&
		  		             (m > r[(i+(2*a))*w+j+(2*b)]) && (m >= r[(i-(2*a))*w+j-(2*b)]) )
		  		        	mid[i*w+j] = 2;	
		  		    }
	  			}
	  		}
	  	midImg.setAllPixels(mid);
	  	JIPImgBitmap resultImg = (JIPImgBitmap)JIPImage.newImage(w, h, image.getType());
	  	resultImg.setAllPixels(r);
	  	return resultImg;
	}
	
	private JIPImgBitmap susanEdgesSmall (JIPImgBitmap image, JIPImgBitmap midImg, int geomThresh, int[] lut) throws JIPException {
		int w=image.getWidth();
		int h=image.getHeight();
		double []pixels=image.getAllPixels();
		double []r=new double[pixels.length];
		double []mid=new double[pixels.length];
		Arrays.fill(mid, 100);
		int n, p, cp;
		double z, c;
		int m, a, b, x, y, s;
		boolean do_symmetry;
		geomThresh = 730; 
		
		for (int i=1;i<h-1;i++)
			for (int j=1;j<w-1;j++) {
				n=100;
		        p=(i-1)*w + j - 1;
		        cp=(int)pixels[i*w+j];

		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-2; 
		        
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=2;
		        n+=lut[256+cp-(int)pixels[p]];
		        p+=w-2; 
		        
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p++]];
		        n+=lut[256+cp-(int)pixels[p]];

				if (n<geomThresh) 
		            r[i*w+j] = geomThresh-n;
		    }

	  	for (int i=2;i<h-2;i++)
	  		for (int j=2;j<w-2;j++) {
	  			if (r[i*w+j]>0)  {
	  		        m=(int)r[i*w+j];
	  		        n=geomThresh - m;
			        cp=(int)pixels[i*w+j];

	  		        if (n>250) {
	  			        p=(i-3)*w + j - 1;
		  		        x=0;y=0;
	
		  		        c=lut[256+cp-(int)pixels[p++]]; x-=c; y-=c;
			            c=lut[256+cp-(int)pixels[p++]]; y-=c;
			            c=lut[256+cp-(int)pixels[p]];   x+=c; y-=c;
			            p+=w-2; 
	
			            c=lut[256+cp-(int)pixels[p]]; x-=c;
			            p+=2;
			            c=lut[256+cp-(int)pixels[p]]; x+=c;
			            p+=w-2;
	
			            c=lut[256+cp-(int)pixels[p++]]; x-=c; y+=c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=c;
			            c=lut[256+cp-(int)pixels[p]]; x+=c; y+=c;
	
		  		        z = Math.sqrt(x*x + y*y);
		  		        if (z > (0.4*n)) /* 0.5 */ {
		  		            do_symmetry=false;
		  		            if (x==0)
		  		            	z=1000000.0;
		  		            else
		  		            	z=(double)y / x;
		  		            if (z < 0) { z=-z; s=-1; }
		  		            else s=1;
		  		            if (z < 0.5) { /* vert_edge */ a=0; b=1; }
		  		            else { if (z > 2.0) { /* hor_edge */ a=1; b=0; }
		  		            else { /* diag_edge */ if (s>0) { a=1; b=1; }
		  		                                   else { a=-1; b=1; }}}
		  		            if ( (m > r[(i+a)*w+j+b]) && (m >= r[(i-a)*w+j-b])) 
		  		            	mid[i*w+j] = 1;
		  		        }
		  		        else
		  		            do_symmetry=true;
		  		    }
		  		    else 
		  		    	do_symmetry=true;
	
	  		        if (do_symmetry) { 
	  			        p=(i-1)*w + j - 1;
		  		        x=0; y=0; s=0;
	
		  		        /*   |      \
		  		             y  -x-  w
		  		             |        \   */

			            c=lut[256+cp-(int)pixels[p++]]; x+=c; y+=c; s+=c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=c;
			            c=lut[256+cp-(int)pixels[p]];   x+=c; y+=c; s-=c;
			            p+=w-3; 
	
			            c=lut[256+cp-(int)pixels[p]];   x+=c;
			            p+=2;
			            c=lut[256+cp-(int)pixels[p]];   x+=c;
			            p+=w-2;
	
			            c=lut[256+cp-(int)pixels[p++]]; x+=c; y+=c; s-=c;
			            c=lut[256+cp-(int)pixels[p++]]; y+=c;
			            c=lut[256+cp-(int)pixels[p]];   x+=c; y+=c; s+=c;
		  		        
		  		        if (y==0)
		  		        	z = 1000000.0;
		  		        else
		  		        	z = (double)x / y;
		  		        if (z < 0.5) { /* vertical */ a=0; b=1; }
		  		        else { 
		  		        	if (z > 2.0) { /* horizontal */ a=1; b=0; }
		  		        	else { /* diagonal */ 
		  		        		if (s>0) { a=-1; b=1; }
		  		        		else { a=1; b=1; }
		  		        	}
		  		        }
		  		        if ( (m > r[(i+a)*w+j+b]) && (m >= r[(i-a)*w+j-b]) )
		  		        	mid[i*w+j] = 2;	
		  		    }
	  			}
	  		}
	  	midImg.setAllPixels(mid);
	  	JIPImgBitmap resultImg = (JIPImgBitmap)JIPImage.newImage(w, h, image.getType());
	  	resultImg.setAllPixels(r);
	  	return resultImg;
	}

	private JIPGeomEdges extractEdges (JIPImgBitmap midImg, int lengthThresh) throws JIPException {
		double[] mid = midImg.getAllPixels();
		int w = midImg.getWidth();
		int h = midImg.getHeight();
		int []p = new int[2];
		
		JIPGeomEdges result = (JIPGeomEdges)JIPImage.newImage(w, h, ImageType.EDGES);
		for (int i=5;i<h-5;i++)
	  		for (int j=5;j<w-5;j++) {
	  			if (mid[i*w+j]<8)  {
	  				mid[i*w+j]=100;
	  				Edge e = new Edge();
  		  			e.addPoint(new Point2D(j,i));
	  				p[0] = i;
	  				p[1] = j;
	  				while (nextPoint(w, h, p, mid)) {
	  					mid[p[0]*w+p[1]] = 100;
	  		  			e.addPoint(new Point2D(p[1],p[0]));
	  				}
	  				e.reverse();
	  				while (nextPoint(w, h, p, mid)) {
	  					mid[p[0]*w+p[1]] = 100;
	  		  			e.addPoint(new Point2D(p[1],p[0]));
	  				}
	  				if (e.length() > lengthThresh)
	  					result.addEdge(e);  
	  			}
	  		}
		
		return result;
	}
	
	/**
	 *Calculates the next point to continue with the trace.<BR>
	 *<ul><B>Input parameters:</B><BR>
	 *<li>rp: Start row<BR>
	 *<li>cp: Start column<BR><BR>
	 *</ul>
	 *<ul><B>Output parameters:</B><BR>
	 *<li>Sets the point in the temporal list ot edges<BR><BR>
	 *</ul>
	 *Additional notes: rp and cp are a class variables, they are used
	 *from TraceEdge to pass some reference values.
	 */
	public boolean nextPoint(int w, int h, int []p, double[]bmp) {
		int i, r, c;
		int[] roff = new int[8];
		int[] coff = new int[8];
		roff[0] = 1;
		roff[1] = 0;
		roff[2] = -1;
		roff[3] = 0;
		roff[4] = 1;
		roff[5] = 1;
		roff[6] = -1;
		roff[7] = -1;
		coff[0] = 0;
		coff[1] = 1;
		coff[2] = 0;
		coff[3] = -1;
		coff[4] = 1;
		coff[5] = -1;
		coff[6] = -1;
		coff[7] = 1;

		for (i = 0; i < 8; i++) {
			r = p[0] + roff[i];
			c = p[1] + coff[i];
			if (r >= 0 && c >= 0 && r < h && c < w)
				if (bmp[c + r * w] < 8) {
					p[0] += roff[i];
					p[1] += coff[i];
					return true;
				}
		}
		return false;
	} 
}